## Measuring Political Support and Issue Ownership Using Endorsement Experiments,
## with Application to Militant Groups in Pakistan

library(R2jags)

J <- 4 ## number of issue areas asked about
K <- 5 ## number of endorsement groups (including control)
n.province <- 4
n.per.province <- c(2247,1280,924,761)
N <- sum(n.per.province) ## number of survey respondents

  ## Ideal Points
  delta.province <- c(-1,-0.25,0.25,1)
  province <- NULL
  for (i in 1:n.province){  province <- c(province, rep(i,n.per.province[i])) }
  
  x <- rep(NA, N)
  for (i in 1:N){ x[i] <- rnorm(1, mean=delta.province[province[i]], sd=1) }
  x <- x / sd(x)
  
  beta <- runif(J, 0.5, 2) ## discrimination parameter
  alpha <- matrix(0.1,J,4)
  for (j in 1:J){  ## variation in thresholds, but distance ensured
	  while(alpha[j,2]-alpha[j,1] < 0.5 | alpha[j,3]-alpha[j,2] < 0.5 | alpha[j,4]-alpha[j,3] < 0.5){
		  tmp <- rnorm(4, mean=0, sd=1)
		  alpha[j,] <- tmp[order(tmp)]
	  }}
  
  theta.province <- matrix(c(rep(0,n.province), c(0.5,0.5,0.5,0.5), c(0,0.5,0.5,1), c(-1,-1,0,1), c(-1,-0.5,-0.5,0)), n.province, K)
  phi <- matrix(c(rep(0,n.province), c(0.3,0.3,0.3,0.3), c(0.3,0.3,0.3,0.3), c(0.3,0.3,0.3,0.3), c(0.3,0.3,0.3,0.3)), n.province, K)

  lambda.province <- array(NA, c(n.province, J, K))
  for (i in 1:n.province){
    for (j in 1:J){
      for (k in 1:K){
        lambda.province[i,j,k] <- rnorm(1, mean=theta.province[i,k], sd=phi[i,k])
  }}}
  
  omega <- matrix(c(rep(0,J), runif(J*(K-1), 0.3, 1)), J, K)
  
  s <- array(NA, c(N,J,K)) ##Individual level impact of group endorsement k
  for (i in 1:N){
    for (j in 1:J){
      for (k in 1:K){
        s[i,j,k] <- rnorm(1, mean=lambda.province[province[i],j,k], sd=omega[j,k])
  }}}

  
##Ordered Responses (5 choices)
Y <- array(NA, c(N,J,K))
for (i in 1:N){
	treat <- ifelse(sample(c(0,1),1)==1,TRUE,FALSE) ## Half of Jake's respondents receive the control, remaining half split among treatments
	for (j in 1:J){
		for (k in 1:K){
			p1 <- pnorm(alpha[j,1] - beta[j]*(x[i] + s[i,j,k]))
			p2 <- pnorm(alpha[j,2] - beta[j]*(x[i] + s[i,j,k])) - (p1)
			p3 <- pnorm(alpha[j,3] - beta[j]*(x[i] + s[i,j,k])) - (p1 + p2)
			p4 <- pnorm(alpha[j,4] - beta[j]*(x[i] + s[i,j,k])) - (p1 + p2 + p3)
			p5 <- 1 - (p1 + p2 + p3 + p4)
			
			p <- c(p1,p2,p3,p4,p5)
			Y[i,j,k] <- sample(c(1:5), size=1, prob=p) ## Full array of potential outcomes
		}}
	if (treat==FALSE) { Y[i,,2:K] <- NA}
	if (treat==TRUE) {
		Y[i,,1] <- NA 
		for (j in 1:J){ 
			answered <- sample(2:K,1); Y[i,j,-answered] <- NA  
		}
	}}

Y2 <- matrix(NA,N,J)
endorser <- matrix(NA, N,J)
for (i in 1:N){
	for (j in 1:J){
		endorser[i,j] <- (1:5)[!is.na(Y[i,j,])]
		Y2[i,j] <- Y[i,j,endorser[i,j]]
	}}
Y <- Y2

  

  lambda.province.real <- lambda.province
  
  #######################################
  ##  Prep and send data through JAGS  ##
  #######################################
  data <- list ("N", "J", "K", "Y", "endorser", "province", "n.province", "lambda.province.real") 
  inits <- function() {list(alpha0=matrix(seq(-2,2,length.out=(4*J)),4,J), beta=runif(J, 0.2, 2), 
    x=rnorm(N), s=matrix(rep(0, N*J), N,J), 
	delta.province=rep(0,n.province),
    lambda.province=array(c(rep(0,n.province*J),rep(0,n.province*J*(K-1))), c(n.province, J, K)), 
    theta.province=matrix(c(rep(0,n.province),rep(0,n.province*(K-1))),n.province,K),
    omega=matrix(c(rep(0,J),runif(J*(K-1))),J,K) )} 
  params <- c("lambda.province.error") 
  model <- "ModelSimulationsIssueProvince.txt" 
  fit <- jags(data, inits, params, n.iter=20000, model.file=model)
  fit <- autojags(fit, n.iter=20000, n.thin=20, Rhat=1.02, n.update=5)
  save.image(file = Outfile)
  
  